#delimit;
clear;
version 12;
set more off, perm;
  quietly log;
  local logon = r(status);
  if "`logon'" == "on" {; log close; };
log using SMPO_MC_SS.log, text replace;

/*******************************/
/*  Mark David Nieman  */
/*  1/21/2016          */
/*  Strategic Logit, Strategic Probit, Split-sample Logit, 
		Split-Sample Probit, Traditional Logit, Traditional Probit  */
/*  DGP: strategic */
/*  Monte Carlos & Comparisons (Clarke distribution-free test) */
/*  Stata 12.1        */
/*******************************/

#delimit;
clear;
clear matrix;
clear mata;
program drop _all;
set seed 1999;


/* Estimators */
	/* strategic logit */
		#delimit;
		program define sl_lf, rclass;
			args lnf A_L A_Rr B_Rr;
			tempvar pB_r pA_R ;
				quietly gen double `pB_r' = (exp(`B_Rr'))/(1+exp(`B_Rr'));
				quietly gen double `pA_R' = (exp(`pB_r'*`A_Rr'))/(exp(`A_L')+exp(`pB_r'*`A_Rr'));			
			quietly replace `lnf' = ln((1-`pA_R')) if $ML_y1==0;
			quietly replace `lnf' = ln((`pA_R')*(1-`pB_r')) if $ML_y1==1;
			quietly replace `lnf' = ln((`pA_R')*(`pB_r')) if $ML_y1==2;
		end;
		
	/* strategic logit with partial observability */
		#delimit;
		program define slpo_lf, rclass;
			args lnf A_L A_Rr B_Rr;
			tempvar pB_r pA_R ;
				quietly gen double `pB_r' = (exp(`B_Rr'))/(1+exp(`B_Rr'));
				quietly gen double `pA_R' = (exp(`pB_r'*`A_Rr'))/(exp(`A_L')+exp(`pB_r'*`A_Rr'));			
			quietly replace `lnf' = ln((1-`pA_R') + (`pA_R')*(1-`pB_r')) if $ML_y1==0;
			quietly replace `lnf' = ln((`pA_R')*(`pB_r')) if $ML_y1==1;
		end;
	
	/* create the split-sample logit */
		#delimit;
		program define ssl_lf, rclass;		
			args lnf gamma beta;
			tempvar rel out;
				quietly gen double `rel' = exp(`gamma')/(1+exp(`gamma'));
				quietly gen double `out' = exp(`beta')/(1+exp(`beta'));
			quietly replace `lnf' = ln((1-`rel')+(`rel'*(1-`out'))) if $ML_y1==0;	
			quietly replace `lnf' = ln((`rel')*(`out')) if $ML_y1==1;
		end;



	#delimit;
/* program to generate y */
program define mc_strat, rclass;
  syntax, indvars(string);
	use "`indvars'.dta", clear;
#delimit;
	/* Agent Error - Gumble Dist (Type 1 Extreme Value) */
	gen double u_r = -ln(-ln(runiform()));
	gen double u_l = -ln(-ln(runiform()));
	gen double u_R = -ln(-ln(runiform()));
	gen double u_L = -ln(-ln(runiform()));
	
	/* Player B Utilities */
	gen double BU_Rl = 0;
	gen double BU_Rr = Z1 + C;
	gen double B_Rl = BU_Rl + u_l;
	gen double B_Rr = BU_Rr + u_r;
	
	/* Expectation of B */
	gen double p_r = exp(B_Rr)/(exp(B_Rl)+exp(B_Rr));
	gen double p_l = exp(B_Rl)/(exp(B_Rl)+exp(B_Rr));
	
	/* Player A Utilities */
	gen double AU_L = X1;
	gen double AU_Rl = 0;
	gen double AU_Rr = X2 + C;
	gen double A_L = AU_L + u_L;
	gen double A_R = p_r*AU_Rr + p_l*AU_Rl + u_R;

	gen double p_R = exp(A_R)/(exp(A_L)+exp(A_R));
	gen double p_L = exp(A_L)/(exp(A_L)+exp(A_R));
	/* Outcomes */
		/* DV for full strategic logit */
		gen     Y3 = . ;
		replace Y3 = 0 if A_L >= A_R;
		replace Y3 = 1 if (A_R > A_L) & (B_Rl >= B_Rr);
		replace Y3 = 2 if (A_R > A_L) & (B_Rr > B_Rl);
		
		/* DV when partially observed */
		gen      Y = . ;
		replace  Y = 0 if (A_L >= A_R) | ((A_R > A_L) & (B_Rl >= B_Rr));
		replace  Y = 1 if (A_R > A_L) & (B_Rr > B_Rl);
		
		/* Summary of each outcome */
		gen out1 = 0;
		replace out1 = 1 if A_L >= A_R;
			sum out1;
			return scalar out1_true=r(mean);
		gen out2 = 0;
		replace out2 = 1 if (A_R > A_L) & (B_Rl >= B_Rr);
			sum out2;
			return scalar out2_true=r(mean);
		gen out3 = 0;
		replace out3 = 1 if (A_R > A_L) & (B_Rr > B_Rl); 
			sum out3;
			return scalar out3_true=r(mean);	
	#delimit;
	/* Compare estimators */
		/* Strategic logit with full information */
		ml model lf sl_lf (A_L: Y3 = X1,nocons) (A_Rr:  = X2 C,nocons) (B_Rr:  = Z1 C,nocons) ;
		ml maximize, diff;	
			return scalar X1_sl_AL = [A_L]_b[X1];
			return scalar X1_sl_AL_se = [A_L]_se[X1];
			return scalar X2_sl_ARr = [A_Rr]_b[X2];
			return scalar X2_sl_ARr_se = [A_Rr]_se[X2];
			return scalar C_sl_ARr = [A_Rr]_b[C];
			return scalar C_sl_ARr_se = [A_Rr]_se[C];
			return scalar Z1_sl_BRr = [B_Rr]_b[Z1];
			return scalar Z1_sl_BRr_se = [B_Rr]_se[Z1];	
			return scalar C_sl_BRr = [B_Rr]_b[C];
			return scalar C_sl_BRr_se = [B_Rr]_se[C];				
			matrix ll_sl = e(ll);
			
			/* Model fit */
				/* i's likelihood */
				gen double pB_r_sl = (exp( [B_Rr]_b[Z1]*Z1 + [B_Rr]_b[C]*C )) 
								  /(1+exp( [B_Rr]_b[Z1]*Z1 + [B_Rr]_b[C]*C ));
				gen double pA_R_sl =     (exp( pB_r_sl*([A_Rr]_b[X2]*X2 + [A_Rr]_b[C]*C) )) 
				/(exp([A_L]_b[X1]*X1 ) +  exp( pB_r_sl*([A_Rr]_b[X2]*X2 + [A_Rr]_b[C]*C) ));
				gen double lnfj_sl = . ;
					replace lnfj_sl = ln((1-pA_R_sl)) if Y3==0;
					replace lnfj_sl = ln((pA_R_sl)*(1-pB_r_sl)) if Y3==1;
					replace lnfj_sl = ln((pA_R_sl)*(pB_r_sl)) if Y3==2;
				
				/* pr(outcome) */
				gen double pr_out1_sl = 1-pA_R_sl;
				gen double pr_out2_sl = pA_R_sl*(1-pB_r_sl);
				gen double pr_out3_sl = pA_R_sl*pB_r_sl;
				gen double predict_out1_sl = 0;
					replace predict_out1_sl = 1 if pr_out1_sl >= max(pr_out3_sl,pr_out2_sl);
					sum predict_out1_sl;
					return scalar predict_out1_sl = r(mean);
				gen double predict_out2_sl = 0;
					replace predict_out2_sl = 1 if pr_out2_sl > max(pr_out1_sl,pr_out3_sl);
					sum predict_out2_sl;
					return scalar predict_out2_sl = r(mean);					
				gen double predict_out3_sl = 0;
					replace predict_out3_sl = 1 if pr_out3_sl > max(pr_out1_sl,pr_out2_sl);
					sum predict_out3_sl;
					return scalar predict_out3_sl = r(mean);					
				/* observations correctly predicted */
				gen correct_sl = 0;
					replace correct_sl = 1 if ( (predict_out1_sl==1 & Y3==0) | (predict_out2_sl==1 & Y3==1) | (predict_out3_sl==1 & Y3==2) );
				sum correct_sl;
					return scalar correct_pct_sl = r(mean);
				/* Y3s correctly predicted */
				gen correct_out3_sl = . ;
					replace correct_out3_sl = 0 if predict_out3_sl==0 & out3==1;
					replace correct_out3_sl = 1 if predict_out3_sl==1 & out3==1;
				sum correct_out3_sl;
					return scalar correct_out3_pct_sl = r(mean);
				/* Y3 false positive */
				gen incorrect_out3_sl = . ;
					replace incorrect_out3_sl = 0 if predict_out3_sl==0 & out3==0;
					replace incorrect_out3_sl = 1 if predict_out3_sl==1 & out3==0;
				sum incorrect_out3_sl;
					return scalar incorrect_out3_pct_sl = r(mean);
				/* Y2s correctly predicted */
				gen correct_out2_sl = . ;
					replace correct_out2_sl = 0 if predict_out2_sl==0 & out2==1;
					replace correct_out2_sl = 1 if predict_out2_sl==1 & out2==1;
				sum correct_out2_sl;
					return scalar correct_out2_pct_sl = r(mean);
				
		#delimit;		
		/* Strategic logit with partial observability */
		ml model lf slpo_lf (A_L: Y = X1,nocons) (A_Rr:  = X2 C,nocons) (B_Rr:  = Z1 C,nocons) ;
		ml maximize, diff;	
			return scalar X1_slpo_AL = [A_L]_b[X1];
			return scalar X1_slpo_AL_se = [A_L]_se[X1];
			return scalar X2_slpo_ARr = [A_Rr]_b[X2];
			return scalar X2_slpo_ARr_se = [A_Rr]_se[X2];
			return scalar C_slpo_ARr = [A_Rr]_b[C];
			return scalar C_slpo_ARr_se = [A_Rr]_se[C];
			return scalar Z1_slpo_BRr = [B_Rr]_b[Z1];
			return scalar Z1_slpo_BRr_se = [B_Rr]_se[Z1];	
			return scalar C_slpo_BRr = [B_Rr]_b[C];
			return scalar C_slpo_BRr_se = [B_Rr]_se[C];				
			matrix ll_slpo = e(ll);
			
			/* Model fit */
				/* i's likelihood */
				gen double pB_r_slpo =   (exp( [B_Rr]_b[Z1]*Z1 + [B_Rr]_b[C]*C )) 
									  /(1+exp( [B_Rr]_b[Z1]*Z1 + [B_Rr]_b[C]*C ));
				gen double pA_R_slpo =   (exp( pB_r_slpo*([A_Rr]_b[X2]*X2 + [A_Rr]_b[C]*C) ))
				  /(exp([A_L]_b[X1]*X1) + exp( pB_r_slpo*([A_Rr]_b[X2]*X2 + [A_Rr]_b[C]*C) ));
				gen double lnfj_slpo = . ;
					replace lnfj_slpo = ln((1-pA_R_slpo) + (pA_R_slpo)*(1-pB_r_slpo)) if Y==0;
					replace lnfj_slpo = ln((pA_R_slpo)*(pB_r_slpo)) if Y==1;
				
				/* pr(outcome) */
				gen double pr_out1_slpo = 1-pA_R_slpo;
				gen double pr_out2_slpo = pA_R_slpo*(1-pB_r_slpo);
				gen double pr_out3_slpo = pA_R_slpo*pB_r_slpo;
				gen double predict_out1_slpo = 0;
					replace predict_out1_slpo = 1 if pr_out1_slpo >= max(pr_out3_slpo,pr_out2_slpo);
					sum predict_out1_slpo;
					return scalar predict_out1_slpo = r(mean);
				gen double predict_out2_slpo = 0;
					replace predict_out2_slpo = 1 if pr_out2_slpo > max(pr_out1_slpo,pr_out3_slpo);
					sum predict_out2_slpo;
					return scalar predict_out2_slpo = r(mean);
				gen double predict_out3_slpo = 0;
					replace predict_out3_slpo = 1 if pr_out3_slpo > max(pr_out1_slpo,pr_out2_slpo);
					sum predict_out3_slpo;
					return scalar predict_out3_slpo = r(mean);
				/* observations correctly predicted (partial obs 0/1) */
				gen correct_slpo = 0;
					replace correct_slpo = 1 if ( (predict_out3_slpo==1 & out3==1) | (predict_out3_slpo==0 & out3==0) );
				sum correct_slpo;
					return scalar correct_pct_slpo = r(mean);					
				/* observations correctly predicted (full info) */
				gen correct_slpo_full = 0;
					replace correct_slpo_full = 1 if ( (predict_out1_slpo==1 & out1==1) | (predict_out2_slpo==1 & out2==1) | (predict_out3_slpo==1 & out3==1) );
				sum correct_slpo_full;
					return scalar correct_pct_slpo_full = r(mean);
				/* Y3s correctly predicted */
				gen correct_out3_slpo = . ;
					replace correct_out3_slpo = 0 if predict_out3_slpo==0 & out3==1;
					replace correct_out3_slpo = 1 if predict_out3_slpo==1 & out3==1;
				sum correct_out3_slpo;
					return scalar correct_out3_pct_slpo = r(mean);
				/* Y3 false positive */
				gen incorrect_out3_slpo = . ;
					replace incorrect_out3_slpo = 0 if predict_out3_slpo==0 & out3==0;
					replace incorrect_out3_slpo = 1 if predict_out3_slpo==1 & out3==0;
				sum incorrect_out3_slpo;
					return scalar incorrect_out3_pct_slpo = r(mean);
				/* Y2s correctly predicted */
				gen correct_out2_slpo = . ;
					replace correct_out2_slpo = 0 if predict_out2_slpo==0 & out2==1;
					replace correct_out2_slpo = 1 if predict_out2_slpo==1 & out2==1;
				sum correct_out2_slpo;
					return scalar correct_out2_pct_slpo = r(mean);		
			
		
		/* Split-population logit */
		gen X1_prime = -X1;
		ml model lf ssl_lf (Rel: Y = X1_prime C,nocons) (Out: Y= X2 Z1 C,nocons);
		ml maximize, diff;
			return scalar X1_ssl_Rel = [Rel]_b[X1_prime];
			return scalar X1_ssl_Rel_se = [Rel]_se[X1_prime];
			return scalar C_ssl_Rel = [Rel]_b[C];
			return scalar C_ssl_Rel_se = [Rel]_se[C];
			return scalar X2_ssl_Out = [Out]_b[X2];
			return scalar X2_ssl_Out_se = [Out]_se[X2];
			return scalar Z1_ssl_Out = [Out]_b[Z1];
			return scalar Z1_ssl_Out_se = [Out]_se[Z1];
			return scalar C_ssl_Out = [Out]_b[C];
			return scalar C_ssl_Out_se = [Out]_se[C];
			matrix ll_ssl = e(ll);
			
			/* Model fit */
				/* i's likelihood */
				quietly gen double rel_ssl =    exp([Rel]_b[X1_prime]*X1_prime + [Rel]_b[C]*C )
											/(1+exp([Rel]_b[X1_prime]*X1_prime + [Rel]_b[C]*C ));
				quietly gen double out_ssl =    exp([Out]_b[X2]*X2 + [Out]_b[Z1]*Z1 + [Out]_b[C]*C )
											/(1+exp([Out]_b[X2]*X2 + [Out]_b[Z1]*Z1 + [Out]_b[C]*C ));
				gen double lnfj_ssl = . ;
					replace lnfj_ssl = ln((1-rel_ssl) + (rel_ssl)*(1-out_ssl)) if Y==0;
					replace lnfj_ssl = ln((rel_ssl)*(out_ssl)) if Y==1;
	
				/* pr(outcome) */
				gen double pr_out1_ssl = 1-rel_ssl;
				gen double pr_out2_ssl = rel_ssl*(1-out_ssl);
				gen double pr_out3_ssl = rel_ssl*out_ssl;
				gen double predict_out1_ssl = 0;
					replace predict_out1_ssl = 1 if pr_out1_ssl >= max(pr_out3_ssl,pr_out2_ssl);
					sum predict_out1_ssl;
					return scalar predict_out1_ssl = r(mean);
				gen double predict_out2_ssl = 0;
					replace predict_out2_ssl = 1 if pr_out2_ssl > max(pr_out1_ssl,pr_out3_ssl);
					sum predict_out2_ssl;
					return scalar predict_out2_ssl = r(mean);
				gen double predict_out3_ssl = 0;
					replace predict_out3_ssl = 1 if pr_out3_ssl > max(pr_out1_ssl,pr_out2_ssl);
					sum predict_out3_ssl;
					return scalar predict_out3_ssl = r(mean);
				/* observations correctly predicted (partial obs 0/1) */
				gen correct_ssl = 0;
					replace correct_ssl = 1 if ( (predict_out3_ssl==1 & out3==1) | (predict_out3_ssl==0 & out3==0) );
				sum correct_ssl;
					return scalar correct_pct_ssl = r(mean);				
				/* observations correctly predicted (full info) */
				gen correct_ssl_full = 0;
					replace correct_ssl_full = 1 if ( (predict_out1_ssl==1 & out1==1) | (predict_out2_ssl==1 & out2==1) | (predict_out3_ssl==1 & out3==1) );
				sum correct_ssl_full;
					return scalar correct_pct_ssl_full = r(mean);
				/* Y3s correctly predicted */
				gen correct_out3_ssl = . ;
					replace correct_out3_ssl = 0 if predict_out3_ssl==0 & out3==1;
					replace correct_out3_ssl = 1 if predict_out3_ssl==1 & out3==1;
				sum correct_out3_ssl;
					return scalar correct_out3_pct_ssl = r(mean);
				/* Y3 false positive */
				gen incorrect_out3_ssl = . ;
					replace incorrect_out3_ssl = 0 if predict_out3_ssl==0 & out3==0;
					replace incorrect_out3_ssl = 1 if predict_out3_ssl==1 & out3==0;
				sum incorrect_out3_ssl;
					return scalar incorrect_out3_pct_ssl = r(mean);
				/* Y2s correctly predicted */
				gen correct_out2_ssl = . ;
					replace correct_out2_ssl = 0 if predict_out2_slpo==0 & out2==1;
					replace correct_out2_ssl = 1 if predict_out2_slpo==1 & out2==1;
				sum correct_out2_ssl;
					return scalar correct_out2_pct_ssl = r(mean);		
											
		/* Traditional logit */
		logit Y X1_prime X2 Z1 C,nocons;
			return scalar X1_l = _b[X1_prime];
			return scalar X1_l_se = _se[X1_prime];
			return scalar X2_l = _b[X2];
			return scalar X2_l_se = _se[X2];
			return scalar Z1_l = _b[Z1];
			return scalar Z1_l_se = _se[Z1];
			return scalar C_l = _b[C];
			return scalar C_l_se = _se[C];
			matrix ll_l = e(ll);
			
			/* Model fit */
				/* i's likelihood */
				gen double lnfj_l = . ;
				replace lnfj_l = ln(invlogit(-(_b[X1_prime]*X1_prime + _b[X2]*X2 + _b[Z1]*Z1 + _b[C]*C ))) if Y==0;
				replace lnfj_l = ln(invlogit(  _b[X1_prime]*X1_prime + _b[X2]*X2 + _b[Z1]*Z1 + _b[C]*C  )) if Y==1;
				
				/* pr(outcome) */
				gen double pr_out3_l = invlogit(_b[X1]*X1 + _b[X2]*X2 + _b[Z1]*Z1 + _b[C]*C );
				gen double predict_out3_l = 0;
					replace predict_out3_l = 1 if pr_out3_l >= .5;
				sum predict_out3_l;
					return scalar predict_out3_l = r(mean);
				/* observations correctly predicted (partial obs 0/1) */
				gen correct_l = 0;
					replace correct_l = 1 if ( (predict_out3_l==1 & out3==1) | (predict_out3_l==0 & out3==0) );
				sum correct_l;
					return scalar correct_pct_l = r(mean);
				/* Y3s correctly predicted */
				gen correct_out3_l = . ;
					replace correct_out3_l = 0 if predict_out3_l==0 & out3==1;
					replace correct_out3_l = 1 if predict_out3_l==1 & out3==1;
				sum correct_out3_l;
					return scalar correct_out3_pct_l = r(mean);
				/* Y3 false positive */
				gen incorrect_out3_l = . ;
					replace incorrect_out3_l = 0 if predict_out3_l==0 & out3==0;
					replace incorrect_out3_l = 1 if predict_out3_l==1 & out3==0;
				sum incorrect_out3_l;
					return scalar incorrect_out3_pct_l = r(mean);

		/* Use Clarke's non-parametric approach to compare each observation's likelihood: slpo vs sl*/
		signtest lnfj_slpo=lnfj_sl;
		return scalar Clark_slpo_sl_pos = r(N_pos);
		return scalar Clark_slpo_sl_neg = r(N_neg);
		return scalar Clark_slpo_sl_tie = r(N_tie);
		return scalar Clark_slpo_sl_pos_1s = r(p_pos);
		return scalar Clark_slpo_sl_neg_1s = r(p_neg);
		return scalar Clark_slpo_sl_eq_2s = r(p_2);
		
		/* Use Clarke's non-parametric approach to compare each observation's likelihood: slpo vs ssl*/
		signtest lnfj_slpo=lnfj_ssl;
		return scalar Clark_slpo_ssl_pos = r(N_pos);
		return scalar Clark_slpo_ssl_neg = r(N_neg);
		return scalar Clark_slpo_ssl_tie = r(N_tie);		
		return scalar Clark_slpo_ssl_pos_1s = r(p_pos);
		return scalar Clark_slpo_ssl_neg_1s = r(p_neg);
		return scalar Clark_slpo_ssl_eq_2s = r(p_2);

		/* Use Clark's non-parametric approach to compare each observation's likelihood: slpo vs logit*/
		signtest lnfj_slpo=lnfj_l;
		return scalar Clark_slpo_l_pos = r(N_pos);
		return scalar Clark_slpo_l_neg = r(N_neg);
		return scalar Clark_slpo_l_tie = r(N_tie);		
		return scalar Clark_slpo_l_pos_1s = r(p_pos);
		return scalar Clark_slpo_l_neg_1s = r(p_neg);
		return scalar Clark_slpo_l_eq_2s = r(p_2);					
		
		/* Vuong Test: slpo vs sl*/
		svmat ll_sl;
		svmat ll_slpo;
		gen num_V_slpo_sl = ll_slpo1 - ll_sl1 + ((5/2)*ln(_N)-(5/2)*ln(_N));
		sum num_V_slpo_sl;
			return scalar num_V_slpo_sl = r(mean);
		egen denom1_V_1_slpo_sl = mean((lnfj_slpo/lnfj_sl)^2);
		egen denom1_V_2_slpo_sl = mean(lnfj_slpo/lnfj_sl);
		gen denom1_V_3_slpo_sl = denom1_V_2_slpo_sl^2;
		gen denom2_V_slpo_sl = sqrt(denom1_V_1_slpo_sl - denom1_V_3_slpo_sl);
		sum denom2_V_slpo_sl;
			return scalar denom2_V_slpo_sl = r(mean);
		gen Vuong_slpo_sl = num_V_slpo_sl/denom2_V_slpo_sl;
		sum Vuong_slpo_sl;
			return scalar Vuong_slpo_sl = r(mean);

		/* Vuong Test: slpo vs ssl*/
		svmat ll_ssl;
		gen num_V_slpo_ssl = ll_slpo1 - ll_ssl1 + ((5/2)*ln(_N)-(5/2)*ln(_N));
		sum num_V_slpo_ssl;
			return scalar num_V_slpo_ssl = r(mean);
		egen denom1_V_1_slpo_ssl = mean((lnfj_slpo/lnfj_ssl)^2);
		egen denom1_V_2_slpo_ssl = mean(lnfj_slpo/lnfj_ssl);
		gen denom1_V_3_slpo_ssl = denom1_V_2_slpo_ssl^2;
		gen denom2_V_slpo_ssl = sqrt(denom1_V_1_slpo_ssl - denom1_V_3_slpo_ssl);
		sum denom2_V_slpo_ssl;
			return scalar denom2_V_slpo_ssl = r(mean);
		gen Vuong_slpo_ssl = num_V_slpo_ssl/denom2_V_slpo_ssl;
		sum Vuong_slpo_ssl;
			return scalar Vuong_slpo_ssl = r(mean);
			
		/* Vuong Test: slpo vs l*/
		svmat ll_l;
		gen num_V_slpo_l = ll_slpo1 - ll_l1 + ((5/2)*ln(_N)-(4/2)*ln(_N));
		sum num_V_slpo_l;
			return scalar num_V_slpo_l=r(mean);
		egen denom1_V_1_slpo_l = mean((lnfj_slpo/lnfj_l)^2);
		egen denom1_V_2_slpo_l = mean(lnfj_slpo/lnfj_l);
		gen denom1_V_3_slpo_l = denom1_V_2_slpo_l^2;
		gen denom2_V_slpo_l = sqrt(denom1_V_1_slpo_l - denom1_V_3_slpo_l);
		sum denom2_V_slpo_l;
			return scalar denom2_V_slpo_l = r(mean);
		gen Vuong_slpo_l = num_V_slpo_l/denom2_V_slpo_l;
		sum Vuong_slpo_l;
			return scalar Vuong_slpo_l = r(mean);			
		
end;

		
	#delimit;
/* Generate Independent Variables */
	set obs 2000;

	/* Generate Covariates U[-4, 4] */
	gen double X1 = runiform()*8 - 4;
	gen double X2 = runiform()*8 - 4;
	gen double X3 = runiform()*8 - 4;
	gen double Z1 = runiform()*8 - 4;
	gen double C  = runiform()*8 - 4;
	
	save indvars.dta, replace;
	
	
/* Run simulations */	
	simulate
		X1_sl_AL=r(X1_sl_AL) X1_sl_AL_se=r(X1_sl_AL_se)
		X2_sl_ARr=r(X2_sl_ARr) X2_sl_ARr_se=r(X2_sl_ARr_se)
		C_sl_ARr=r(C_sl_ARr) C_sl_ARr_se=r(C_sl_ARr_se)
		Z1_sl_BRr=r(Z1_sl_BRr) Z1_sl_BRr_se=r(Z1_sl_BRr_se)
		C_sl_BRr=r(C_sl_BRr) C_sl_BRr_se=r(C_sl_BRr_se)
		
		X1_slpo_AL=r(X1_slpo_AL) X1_slpo_AL_se=r(X1_slpo_AL_se)
		X2_slpo_ARr=r(X2_slpo_ARr) X2_slpo_ARr_se=r(X2_slpo_ARr_se)
		C_slpo_ARr=r(C_slpo_ARr) C_slpo_ARr_se=r(C_slpo_ARr_se)
		Z1_slpo_BRr=r(Z1_slpo_BRr) Z1_slpo_BRr_se=r(Z1_slpo_BRr_se)
		C_slpo_BRr=r(C_slpo_BRr) C_slpo_BRr_se=r(C_slpo_BRr_se)
		
		X1_ssl_Rel=r(X1_ssl_Rel) X1_ssl_Rel_se=r(X1_ssl_Rel_se)
		C_ssl_Rel=r(C_ssl_Rel) C_ssl_Rel_se=r(C_ssl_Rel_se)
		X2_ssl_Out=r(X2_ssl_Out) X2_ssl_Out_se=r(X2_ssl_Out_se)
		Z1_ssl_Out=r(Z1_ssl_Out) Z1_ssl_Out_se=r(Z1_ssl_Out_se)
		C_ssl_Out=r(C_ssl_Out) C_ssl_Out_se=r(C_ssl_Out_se)
		
		X1_l=r(X1_l) X1_l_se=r(X1_l_se)
		X2_l=r(X2_l) X2_l_se=r(X2_l_se)
		Z1_l=r(Z1_l) Z1_l_se=r(Z1_l_se)
		C_l=r(C_l) C_l_se=r(C_l_se)
		
		out1_true=r(out1_true)
		out2_true=r(out2_true)
		out3_true=r(out3_true)
		
		predict_out1_sl = r(predict_out1_sl)
		predict_out2_sl = r(predict_out2_sl)
		predict_out3_sl = r(predict_out3_sl)
		predict_out1_slpo = r(predict_out1_slpo)
		predict_out2_slpo = r(predict_out2_slpo)
		predict_out3_slpo = r(predict_out3_slpo)
		predict_out1_ssl = r(predict_out1_ssl)
		predict_out2_ssl = r(predict_out2_ssl)
		predict_out3_ssl = r(predict_out3_ssl)	
		predict_out3_l = r(predict_out3_l)		
		
		correct_pct_sl=r(correct_pct_sl) 
		correct_out3_pct_sl=r(correct_out3_pct_sl)
		incorrect_out3_pct_sl=r(incorrect_out3_pct_sl)
		correct_out2_pct_sl=r(correct_out2_pct_sl)
		
		correct_pct_slpo=r(correct_pct_slpo)
		correct_pct_slpo_full=r(correct_pct_slpo_full)
		correct_out3_pct_slpo=r(correct_out3_pct_slpo)
		incorrect_out3_pct_slpo=r(incorrect_out3_pct_slpo)
		correct_out2_pct_slpo=r(correct_out2_pct_slpo)
		
		correct_pct_ssl=r(correct_pct_ssl)
		correct_pct_ssl_full=r(correct_pct_ssl_full)
		correct_out3_pct_ssl=r(correct_out3_pct_ssl)
		incorrect_out3_pct_ssl=r(incorrect_out3_pct_ssl)
		correct_out2_pct_ssl=r(correct_out2_pct_ssl)
		
		correct_pct_l=r(correct_pct_l)
		correct_out3_pct_l=r(correct_out3_pct_l)
		incorrect_out3_pct_l=r(incorrect_out3_pct_l)
		
		Clark_slpo_sl_pos=r(Clark_slpo_sl_pos)
		Clark_slpo_sl_pos_1s=r(Clark_slpo_sl_pos_1s)
		Clark_slpo_sl_neg_1s=r(Clark_slpo_sl_neg_1s)
		Clark_slpo_sl_eq_2s=r(Clark_slpo_sl_eq_2s)
		
		Clark_slpo_ssl_pos=r(Clark_slpo_ssl_pos)
		Clark_slpo_ssl_pos_1s=r(Clark_slpo_ssl_pos_1s)
		Clark_slpo_ssl_neg_1s=r(Clark_slpo_ssl_neg_1s)
		Clark_slpo_ssl_eq_2s=r(Clark_slpo_ssl_eq_2s)
		
		Clark_slpo_l_pos=r(Clark_slpo_l_pos)
		Clark_slpo_l_pos_1s=r(Clark_slpo_l_pos_1s)
		Clark_slpo_l_neg_1s=r(Clark_slpo_l_neg_1s)
		Clark_slpo_l_eq_2s=r(Clark_slpo_l_eq_2s)
		
		num_V_slpo_sl=r(num_V_slpo_sl)
		denom2_V_slpo_sl=r(denom2_V_slpo_sl)
		Vuong_slpo_sl=r(Vuong_slpo_sl)
		
		num_V_slpo_ssl=r(num_V_slpo_ssl)
		denom2_V_slpo_ssl=r(denom2_V_slpo_ssl)
		Vuong_slpo_ssl=r(Vuong_slpo_ssl)
		
		num_V_slpo_l=r(num_V_slpo_l)
		denom2_V_slpo_l=r(denom2_V_slpo_l)
		Vuong_slpo_l=r(Vuong_slpo_l)		
	,
	
	reps(1000): 
	mc_strat, indvars(indvars) ;
	
	save mc_ss_results, replace;
	
/* View Results */
sum;

	 /*generate root mean standard error for each parameter */
	#delimit;
	gen rmse_X1_sl_AL = sqrt((X1_sl_AL - 1)^2 + (X1_sl_AL_se));
	gen rmse_X1_slpo_AL = sqrt((X1_slpo_AL - 1)^2 + (X1_slpo_AL_se));
	gen rmse_X1_ssl_Rel = sqrt((X1_ssl_Rel - 1)^2 + (X1_ssl_Rel_se));
	gen rmse_X1_l = sqrt((X1_l - 1)^2 + (X1_l_se));
	
	
	gen rmse_X2_sl_ARr = sqrt((X2_sl_ARr - 1)^2 + (X2_sl_ARr_se));
	gen rmse_X2_slpo_ARr = sqrt((X2_slpo_ARr - 1)^2 + (X2_slpo_ARr_se));
	gen rmse_X2_ssl_Out = sqrt((X2_ssl_Out - 1)^2 + (X2_ssl_Out_se));
	gen rmse_X2_l = sqrt((X2_l - 1)^2 + (X2_l_se));

	gen rmse_C_sl_ARr = sqrt((C_sl_ARr - 1)^2 + (C_sl_ARr_se));
	gen rmse_C_slpo_ARr = sqrt((C_slpo_ARr - 1)^2 + (C_slpo_ARr_se));
	gen rmse_C_ssl_Rel = sqrt((C_ssl_Rel - 1)^2 + (C_ssl_Rel_se));

	gen rmse_Z1_sl_BRr = sqrt((Z1_sl_BRr - 1)^2 + (Z1_sl_BRr_se));
	gen rmse_Z1_slpo_BRr = sqrt((Z1_slpo_BRr - 1)^2 + (Z1_slpo_BRr_se));
	gen rmse_Z1_ssl_Out = sqrt((Z1_ssl_Out - 1)^2 + (Z1_ssl_Out_se));
	gen rmse_Z1_l = sqrt((Z1_l - 1)^2 + (Z1_l_se));

	gen rmse_C_sl_BRr = sqrt((C_sl_BRr - 1)^2 + (C_sl_BRr_se));
	gen rmse_C_slpo_BRr = sqrt((C_slpo_BRr - 1)^2 + (C_slpo_BRr_se));
	gen rmse_C_ssl_Out = sqrt((C_ssl_Out - 1)^2 + (C_ssl_Out_se));
	gen rmse_C_l = sqrt((C_l - 1)^2 + (C_l_se));

sum rmse*;
	
/* graph the betas */
#delimit;
twoway kdensity X1_sl_AL, scheme(s1color) color(black) lwidth(thin) lpattern(solid)
	|| kdensity X1_slpo_AL, color(black) lwidth(thick) lpattern(solid)
	|| kdensity X1_ssl_Rel, color(gs5) lwidth(medium) lpattern(dash)
	|| kdensity X1_l, color(gs9) lwidth(medium) lpattern(longdash_dot)
	xtitle("") xline(1, lcolor(gs8) lpattern(dash)) title("X1 (AL)")
	legend(label(1 "Full Info. Strat. Logit") label(2 "Strat. Logit Part. Obs.") 
		   label(3 "Split-Sample Logit") label(4 "Traditional Logit") symx(*.5))
	ylabel("");
	graph save Strat_x1_AL_compare_ss, replace;
	

#delimit;
twoway kdensity X2_sl_ARr, scheme(s1color) color(black) lwidth(thin) lpattern(solid)
	|| kdensity X2_slpo_ARr, color(black) lwidth(thick) lpattern(solid)
	|| kdensity X2_ssl_Out, color(gs5) lwidth(medium) lpattern(dash)
	|| kdensity X2_l, color(gs9) lwidth(medium) lpattern(longdash_dot)
	xtitle("") xline(1, lcolor(gs8) lpattern(dash)) title("X2 (ARr)")
	legend(label(1 "Full Info. Strat. Logit") label(2 "Strat. Logit Part. Obs.") 
		   label(3 "Split-Sample Logit") label(4 "Traditional Logit") symx(*.5))
	ylabel("");
	graph save Strat_x2_ARr_compare_ss, replace;

#delimit;
twoway kdensity C_sl_ARr, scheme(s1color) color(black) lwidth(thin) lpattern(solid)
	|| kdensity C_slpo_ARr, color(black) lwidth(thick) lpattern(solid)
	|| kdensity C_ssl_Rel, color(gs5) lwidth(medium) lpattern(dash)
	xtitle("") xline(1, lcolor(gs8) lpattern(dash)) title("Xc (ARr)")
	legend(label(1 "Full Info. Strat. Logit") label(2 "Strat. Logit Part. Obs.") 
		   label(3 "Split-Sample Logit")  symx(*.5))
	ylabel("");
	graph save Strat_c_ARr_compare_ss, replace;
	
#delimit;
twoway kdensity Z1_sl_BRr, scheme(s1color) color(black) lwidth(thin) lpattern(solid)
	|| kdensity Z1_slpo_BRr, color(black) lwidth(thick) lpattern(solid)
	|| kdensity Z1_ssl_Out, color(gs5) lwidth(medium) lpattern(dash)
	|| kdensity Z1_l, color(gs9) lwidth(medium) lpattern(longdash_dot)
	xtitle("") xline(1, lcolor(gs8) lpattern(dash)) title("X3 (Br)")
	legend(label(1 "Full Info. Strat. Logit") label(2 "Strat. Logit Part. Obs.") 
		   label(3 "Split-Sample Logit") label(4 "Traditional Logit") symx(*.5))
	ylabel("");
	graph save Strat_z1_BRr_compare_ss, replace;
	
#delimit;
twoway kdensity C_sl_BRr, scheme(s1color) color(black) lwidth(thin) lpattern(solid)
	|| kdensity C_slpo_BRr, color(black) lwidth(thick) lpattern(solid)
	|| kdensity C_ssl_Out, color(gs5) lwidth(medium) lpattern(dash)
	|| kdensity C_l, color(gs9) lwidth(medium) lpattern(longdash_dot)
	xtitle("") xline(1, lcolor(gs8) lpattern(dash)) title("Xc (Br)")
	legend(label(1 "Full Info. Strat. Logit") label(2 "Strat. Logit Part. Obs.") 
		   label(3 "Split-Sample Logit") label(4 "Traditional Logit") symx(*.5))
	ylabel("");
	graph save Strat_c_BRr_compare_ss, replace;	
	
	
#delimit;
grc1leg Strat_x1_AL_compare_ss.gph Strat_x2_ARr_compare_ss.gph Strat_c_ARr_compare_ss.gph
		Strat_z1_BRr_compare_ss.gph Strat_c_BRr_compare_ss.gph,
	scheme(s1color) xcommon col(3) holes(4)
	title("")	l1title("Kernel Density")
	leg(Strat_x1_AL_compare_ss.gph)
	b1title("Coefficient") imargin(tiny)
	note("Note: Dashed vertical line represents the equation's true coefficient. Results of 1000 simulations with" 
		 "2000 observations each. Bij represents the estimated coeffient for regressor Xij, where i is the actor"
		 "and j is the outcome. Xc is a common regressor that appears in each actor's utility. Because the"
		 "traditional and split-sample logit models estimate less parameters, they do not appear in each graph.");
		graph save Strat_Coef_compare_ss, replace;
		graph export Strat_Coef_compare_ss.eps, replace;	

		log off;
